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(N ■ 

, ^ ■ We demonstrate a possibility to generate localized states in effectively one-dimensional Bose- 

^ ' Einstein condensates with a negative scattering length in the form of a dark soliton in the presence 

, of an optical lattice (OL) and/or a parabolic magnetic trap. We connect such structures with twisted 

^ ■ localized modes (TLMs) that were previously found in the discrete nonlinear Schrodinger equation. 

' Families of these structures are found as functions of the OL strength, tightness of the magnetic 

CO , trap, and chemical potential, and their stability regions are identified. Stable bound states of two 

TLMs are also found. In the case when the TLMs are unstable, their evolution is investigated by 
means of direct simulations, demonstrating that they transform into large-amplitude fundamental 
solitons. An analytical approach is also developed, showing that two or several fundamental solitons, 
^ ' with the phase shift tt between adjacent ones, may form stable bound states, with parameters quite 

close to those of the TLMs revealed by simulations. TLM structures are found numerically and 

■ explained analytically also in the case when the OL is absent, the condensate being confined only 
by the magnetic trap. 

O ' I. INTRODUCTION 

The current experimental and theoretical studies of Bose-Einstein condensates (BECs) [1] have attracted a great deal 
of interest to nonlinear patterns which can exist in them, including dark [2] and bright [3] solitons. Two-dimensional 
^ (2D) excitations in the form of vortices were also realized experimentally [4]. Other nonlinear excitations, such as 
Faraday waves [5], ring dark solitons and vortex necklaces [6] were predicted to occur in BECs. 
\^ ' The behavior of the condensate crucially depends on the sign of the atomic interactions: dark (bright) solitons can 
. be created in BECs with repulsive (attractive) interactions, resulting from the positive (negative) scattering length. 
' Recent experiments have demonstrated that "tuning" of the effective scattering length, including a possibility to 

■ change its sign, can be achieved by means of a Feshbach resonance [7]. 

' In this work we demonstrate that another type of localized nonlinear excitations, which may be regarded as dark 
"j^ solitons embedded in bright ones, can be created in attractive BECs. They may also be considered as bound states 
of two bright solitons with a phase shift tt between them. In some respects, we follow the path of very recent results 
I of Ref. [8]. These, in turn, were based on earlier works on the so-called twisted localized modes (TLMs) in discrete 
systems [9]; that is why we apply the same term, TLM, to these solitons of the combined dark-bright type. TLMs in 
5 . a discrete system are bright solitons in which the field vanishes, changing its sign, at the central point. An alternative 
way to view a TLM is as a concatenation of an "up-pulse" (i.e., sech(a; — xi)) with a "down pulse" (i.e., — sech(a; — X2)), 
where xi and X2 are appropriately separated, to form an up-down, two-pulse configuration. Notice, however, that our 
work is essentially different from Ref. [8] , since we focus on the case of BECs with negative scattering length (such as 
^Li [10] and ®^Rb [11]), and demonstrate robustness and stability of TLMs in the latter context. These structures, if 
?H regarded as effectively dark solitons in the attractive BEC, are the inverse of recently predicted [12] bright solitons in 
repulsive BECs with an optical lattice. It is expected that a characteristic longitudinal length of the TLM soliton in 
physical units will be on the order of the wavelength of light which induces the OL, i.e., ~ 1 fim, and the number of 
atoms in the soliton may typically be ^ 10'^ - 10^. 

We consider a quasi-lD (cigar-shaped) BEC with negative scattering length, the corresponding normalized Gross- 
Pitaevskii (GP) equation being [1,13,14] 

iut = —Uxx — \u\^u -I- \ex'^ ^' Vo cos^ (27ra;/A)] w, (1) 

where u{x, t) is the macroscopic wave function, e measures the strength of the parabolic magnetic trap, and the last 
term accounts for the OL potential created by interference of two counterpropagating coherent light beams [15-20]. 
As we consider only (effectively) ID configurations, the attractive interactions between the atoms (characterized by 
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the negative scattering length) can not produce collapse in the framework of the present model, hence various solitons 
have a chance to be stable. 

It is obvious that scaling invariance of Eq. (1) makes it possible to fix the optical wavelength, setting k = 27r/A = 1. 

This essentially means that all the remaining length scales of the problem arc measured by comparison to the period 
of the optical lattice which has been fixed to tt. We follow this convention below. The number of atoms in the 
condensate, given by 

/ + 00 
\u{x)\^dx, (2) 
-oo 

is a dynamical invariant of Eq. (1). 

In the next section we report detailed numerical results for TLM solutions in attractive BECs, as modeled by Eq. 
(1). In section III, we will support some of the numerical findings by analytical considerations. Section IV concludes 
the paper. 



II. NUMERICAL RESULTS 



A. Stationary twisted-localized-mode solutions and their linear stability 

We seek stationary solutions to Eq. (1) in the form, 

u{x,t) = ex.p{—i^t)v{x), 
where is the chemical potential, and a real function v{x) obeys the equation 

jjv = —Vxx — + [ex^ + Vo cos^(27rx/A)] w, (3) 

supplemented with free boundary conditions. We have checked that the boundary conditions do not significantly afi'ect 
the results. The corresponding boundary-value problem was solved by means of a finite-difference discretization. A 
second order difference scheme with spacing Ax = 0.2 was used to approximate Vxx- The resulting set of nonlinear 
algebraic equations was solved by means of a Newton iteration. More details on finite-difference schemes and Newton 
type methods can be found in [21]. The solutions of interest were obtained by choosing, as the initial condition for 
the Newton iterations, the following functional form: 

<x) = f^(-l)^7?sech (^-^{x - ,v^V^. (4) 



The ansatz (4) approximates a superposition of nonlinear-Schrodinger (NLS) solitons. The presence of the OL suggests 
choosing positions of the solitons' centers, ^j , as = j7r/2, with some integer j. The sign factor (—1)-' implies the 
phase difference tt between adjacent solitons. This choice of the phase pattern is suggested by the known result that, 
in the discrete NLS equation, only such a configuration may be stable, see Ref. [22] and references therein. Similar 
results for a BEC model with elliptic- function potentials were obtained in Ref. [23]. The number n of solitons in the 
initial ansatz (4) was taken as n = 2, in order to find the TLM soliton proper, or n = 3, with the objective to find a 
bound state of two TLMs, see below. However, it will become clear below that concatenating more of such "building 
block" structures (such as the ones with n = 2 or n = 3) one can, in principle, construct multi-pulse configurations 
with a larger number (n) of elementary pulses. 

Once TLM states were found as solutions to Eq. (3) [starting with the initial approximation (4)], their stability 
was subsequently analyzed by examining the corresponding linearized problem for small perturbations. To this end, 
the perturbed solution was taken as 

u{x, t) = e-'^* ^^v{x) + 5 [a{x)e-^* + 6(a;)e''"**] } , (5) 

where S is an infinitesimal amplitude of the perturbation, a{x) and b{x) are eigenfunctions of a perturbation mode, 
and U! is the corresponding eigenfrequency; * standing for the complex conjugation. Using the ansatz of Eq. (5) to 
0{6), one obtains the linearization problem as an eigenvalue problem where oj is the corresponding eigenfrequency, 
while {a,b*y (where the T denotes transpose), is the corresponding eigenfunction. Using the same finite difference 
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scheme as the one used for computing v(x), we convert this to a matrix eigenvahic problem and find all the eigenvahies 
of the corresponding matrix. This is done by standard matrix eigenvalue solvers built into Matlab [24]. lu can, in 
general, be complex, i.e., u = Wr+itJi, where the subscripts denote the real and imaginary part of the eigenfrequency. 
If cigcnfrcqucncies with non-zero imaginary parts exist, they lead to exponential instabilities, while their absence 
implies linear stability of the solution under consideration. This is monitored by the spectral plane pictures of the 
imaginary (ui) versus the real {u)r) part of the eigenfrequency (see e.g., the bottom right panels of Figs. 1-3 below). 
The presence of even a single lo with LUi ^ implies instability. In all the cases when stable stationary solitons were 
identified according to this criterion (see below), the stability was also checked and confirmed by direct simulations 
of the full equation (1). For more details on the setup and solution of the eigenvalue problem, the interested reader 
is referred to [25,26]. 

To systematically trace the evolution of numerically exact TLM solutions of Eq. (1), we performed one-parameter 
continuations in the three-dimensional parameter space of Eq. (3), e, Vq). Recall that we have set A = 27r, hence 
A is not a free parameter. We start by varying the OL-potential's strength Vq for fixed e = and /i = —1.5. We note 
that negative values of the chemical potential are typical of cases when the corresponding wave-function pattern is 
localized, although ^ can sometimes be positive (see below), which is explained by the presence of the OL potential 
in Eq. (3). 

Numerical results corresponding to this case (varying Vq) are presented in Fig. 1. In the top subplot of the upper 
part of the figure, the number of atoms in the soliton, defined as per Eq. (2), is shown versus Vq, and the largest 
instability growth rate, found from the associated linearized problem, is shown in the bottom panel of the upper part 
of Fig. 1 (zero instability growth rate implies that the soliton is stable). Three characteristic examples of the TLM 
solution are displayed, together with their linear-stability eigenfrequencies, in the lower part of Fig. 1. 

This solution branch exists only for Vq > 0.075. The presence of this cutoff point is expected, as the optical lattice 
is crucial for the existence of TLM solutions, which are not found in the usual NLS equation (at least, in the absence 
of the parabolic-trap potential). Nevertheless, below we will present a case in which a TLM/dark soliton solution is 
possible even in the absence of OL, being supported solely by the magnetic trap. In the interval 0.275 < Vq < 0.7, this 
soliton family is subject to oscillatory instability which is accounted for by the so-called Hamiltonian Hopf bifurcation 
[27]. In the latter, the collision of two pairs of eigenvalues of opposite Krein signature [27] results in the generation of 
a quartet of genuinely complex eigenvalues. The strongest instability, with Imw w 0.107 occurs at Vq « 0.5. 
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FIG. 1. Upper part: the top panel shows the number of particles vs. the OL-potential's strength. The bottom panel shows 
the largest imaginary part of unstable eigenfrequencies (i.e., the instability growth rate) as a function of Vq. Lower part: 
examples of the TLM solutions for Vb = 1 (top panel), Vb = 0.5 (middle panel) and Vo — 0.1 (bottom panel) are shown by 
dashed curves in the left parts of the panels, while the solid line represents the OL potential. The right panel shows the spectral 
plane, uirjiOi), of the real and imaginary parts of the corresponding eigenfrequencies. 



Similarly, we examined the evolution of the solutions, setting Vb = 1 and = — 1.5 in Eq. (3) and performing 
continuation in e, starting from e = 0. In this case too, oscillatory instability was found in a finite interval, 0.05 < 
e < 1.61. The maximum instability growth rate (which is very large in this case) is 0.48, and occurs at e w 0.6. The 
results for this case are shown in Fig. 2. 
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FIG. 2. The same as Fig. 1, but versus the variation of the parabolic-trap strength e. In the lower part, three examples of 
the solutions and their stability are given for e = 0, e = 0.6 and e = 1.75. 



In these cases, the oscillatory instability arises quite naturally. Indeed, similar instabilities had been observed in 
the discrete counterpart of the system and, in fact, in the next section we will give a theoretical justification for their 
occurrence. Another interesting numerical finding is that the TLM solutions appear to be much more robust than 
one would initially suspect. In particular, there seem to exist TLM solutions even when the OL has been almost 
completely smeared out by the harmonic trap (see e.g., the last panel of Fig. 2). We will return to this point in the 
following section. 

We also performed the continuation of the TLM solutions with respect to the chemical potential /j., which has 
yielded results displayed in Fig. 3. In this case, e = and Vq = 1 were chosen. The solutions exist in the region 
M < /^max ~ 0.125. The fact that the cutoff value /imax is positive is explained by a contribution of the mean 
value of the OL potential VqCOS^x in Eq. (3). In this case, oscillatory instability is found in two finite intervals, 
-0.375 <ii< -O.I, and -0.075 < ;U < 0. 
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FIG. 3. The same as in Fig. 1 for tlie variation witli respect to tlie chemical potential ji. In the lower part, three examples 
of the TLM solutions correspond to ji = —10 (top), ji = —0.25 (middle) and = 0.125 (bottom). 



B. Bound states of two TLM solitons 



Complexes which may be considered as a bound state of two TLMs, or simply as a concatenation of three funda- 
mental solitons with phase difference tt between adjacent onc's. wcirc; also studied. In this case, for symmetry purposes, 
cos a; in Eqs. (1) and (3) was replaced by sinx (recall we have set A = 27r), and three pulses in the initial ansatz (4) 
were set at = — tt, 0, tt. We performed the continuation in e, which has demonstrated that the bound state persists 
to very large values of e, see Fig. 4. The branch becomes unstable, through a quartet of complex eigenfrequencies, 
at e > 0.015. A second instability, accounted for by another eigenvalue quartet, sets in at e > 0.075. This second 
oscillatory instability can be easily explained: as we argue below, for each pair of anti-phase pulses there exists a pair 
of eigenvalues with negative Krein signature, which is known to give rise to a Hamiltonian Hopf bifurcation [27]. The 
strongest instability is found at e = 0.23, when the two instability growth rates are 0.39 and 0.35. The first quartet 
returns to stability at e > 0.29, while the second quartet is stabilized for e > 0.33. The bound states are completely 
stable for e > 0.33. 
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FIG. 4. The same as in Fig. 2, but for the bound state of two TLM solitons. In the lower part, three examples are shown, 
for e = (top), e = 0.23 (middle) and e = 2 (bottom panel). Notice that in the second panel of the figure two sets of unstable 
eigenfrequencies are shown (one by stars and one by circles). This indicates the presence of two unstable eigenvalue quartets, 
as can also be verified by the spectral plane figure corresponding to the case of e = 0.23; see also the discussion in the text. 



C. Evolution of unstable TLM solitons 



In the cases for which the TLM solitons or their bound states were found to be unstable (due to the Hamiltonian 
Hopf bifurcations), we simulated nonlinear development of the instability in the framework of the full GP equation 
(1). To this end, a finite-difference scheme with dx = 0.2 [the same stepsize as in dealing with the stationary equation 
(3)] was used, and time integration was performed by means of the fourth-order Rungc-Kutta scheme with a time 
step dt = 0.001. Three examples of the evolution of unstable configurations are shown in Fig. 5. The two top ones 
correspond, respectively, to the strong and weak oscillatory instability of the TLM soliton. In both of these cases, the 
TLM state; evolves into a single fundamental soliton. A difference between the two cases is that, in the former one, 
with the relatively weak OL, the newly formed soliton keeps moving between two minima of the potential, while in 
the latter case , which corresponds to a stronger OL, the fundamental soliton becomes pinned. Finally, the bottom 
subplot of Fig. 5 demonstrates that the oscillatory instability of a bound state of two TLM solitons transforms it into 
a single fundamental soliton with a large amplitude. 
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FIG. 5. Typical examples of the directly simulated evolution of unstable TLM-type solutions. In each set, the top panel 
displays evolution of the negative local density, —\u{x,t)\^ (using a gray-scale space-time plot), the middle panel shows the time 
evolution of the maxiunim of \u{x,t)\^, and the spatial distribution of the density at the end of the simulation is presented in 
the bottom panel. In the latter panel, |w(x,t)|^ in final and initial states, and the potential are shown, respectively, by solid, 
dashed, and dashed-dotted lines. The first example is displayed for e = 0.6 and Vo = 1; recall that the instability growth rate 
is ~ 0.48 in this case, i.e., it is a case of strong instability. A random-noise perturbation of an amplitude 10"' was added to 
the initial configuration in order to excite the instability. The second example corresponds to e = and Vo = 1 for /i — —0.25. 
As the instability growth rate has a very small value in this case, 0.032, a random-perturbation seed with a larger amplitude. 
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was used. The last example pertains to an unstable bound state of two TLM solitons in the case of e = 0.23 and 



Vo = 1.0. In this case, a uniformly distributed random perturbation of amplitude 10 was used to excite the instability. 
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III. ANALYTICAL RESULTS 



A. The minimum size of the TLM state 



The existence of the TLM state in the GP equation with OL can be explained, by means of perturbation theory, as 
a bound state of two fundamental solitons with phase difference tt between them. Here, we consider the case of e = 
(no parabolic trap), but the analysis can be easily generalized to include the trap. If the present model is considered 
as a perturbed NLS equation, we can use an ansatz in the form of Eq. (4) for two solitons. The system's Hamiltonian 
reads 



H 



f 

J — < 



1 



ul^ + Focos^ 



27ra;\ , ,9 



dx. 



Recall that we have set e = 0; here, for generality's sake, we do not fix A = 27r. Then, known methods (see a review 
in Ref. [28]) yield an effective net potential of the interaction between two solitons and interaction of both solitons 
with the OL. This potential is realized as a part of the Hamiltonian that depends on the coordinates ^1 and ^2 of the 
solitons and phase difference At^ between them (cf. a similar potential for solitons in the discrete-NLS model found 
in Ref. [22]) and is of the form: 



Kff (6,6; ^4') = -8?7^ cos (A(/)) exp 



V2 



+- 



2V27r2Fo 



Asinh 



/2\/22f.A 
\ Xr, ) I 



/ 4^277. 

cos I ^ 4l 



cos ^6 



(6) 



The second term in this expression is the Peierls-Nabarro potential (see, e.g., Ref. [25] and references therein) induced 
by the OL. Equilibrium positions of the two-soliton system are determined by equations 



0^2 5(A<^) 



0. 



(7) 



Straightforward consideration of Eqs. (7) and (6) shows that a bound state of two solitons with A0 = tt, which 
corresponds to a TLM-like pattern, may exist if the separation L = — between them exceeds the minimum 
value. 



V2 
V 



In 



2^4 



V27r3Vo 



sinh 



'2^/2^ 
Xt] 



(8) 



This prediction for the minimum separation between fundamental-soliton components of the TLM state was tested 
numerically. As a typical example, in Fig. 6 we show the variation of Z/min with Vb for A = 27r and 77 = \/3. The 
distance between the centers of the solitons was extracted from numerical data as the distance between two local 
maxima of |w(a;)|, an inherent discretization error of ±0.2 being imposed by the finite-difference scheme. As is seen, 
the agreement between the theoretical prediction and the numerically approximated values of Z/min is quite good. 

A more detailed analytical consideration of dynamical properties of a perturbed bound state (following the lines 
of Ref. [22]) shows that, within the framework of the present approximation, the bound state is stable. Moreover, 
the existence of stable bound states of three fundamental solitons, which correspond to the bound state of two TLMs 
considered above, can also be demonstrated by means of this approach. 



9 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 

FIG. 6. The analytical prediction (8) for the niinirnuni separation between two 7r-out-of-phasc fundamental solitons in the 
bound state (stars), and the separation between two density maxima in the TLM found from numerical data (circles), versus 
Vo- Notice that the numerical results should be considered with an error bar of ±0.2, due to the discretization used. 



B. Oscillatory instabilities of the single-TLM and double-TLM states 

It is well known (see e.g., Ref. [25] and references therein) that the Hnearized equations for the perturbations defined 
in Eq. (5) can be rewritten, in terms of new variables U = a + b* and W = a — b*, as 

ujU = L_W=-W^:,-v'^W + {V{x)-h)W, (9) 
ujW = L+U = -U^^ - 3v^U + {V{x) - fx)Un. (10) 

The Krein signature of the eigenvalue can be written as [26] 

K = sign ^- J {UW* + WU*)dx^ . 

Analytical consideration is possible for w's that are real; then U and W are also real, hence K amounts to the 
expression 

K = sign(^-J UWdx^ . (11) 

As was done above, we approximate a TLM solution as a superposition, v{x) = P{x) + N{x), of two stationary 
fundamental-soliton waveforms, positive P{x) and negative N{x), the underlying assumption being a weak overlap 
between P{x) and N{x). Since P and N are, essentially, individual solitons, in the lowest-order approximation they 
obey equations [V{x) - ^] P = P^^ + p3 and [V{x) - fj]N = N^^ + N^. Summing them up, and taking into account 
that the product |PA''| is everywhere much smaller than P^ + (due to the assumed weak overlap between P 
and N), we conclude that, up to the same lowest-order accuracy, the function P(x) + N{x) satisfies the equation 
L-(P + N) = [recall that the operator L_ is defined in Eq. (9)]. Thus, in the lowest approximation, P(x) + N{x) is 
a zero mode of the operator L_. In a similar way, the following approximate result can be obtained for the function 
P{x) - N{xy. 

L-{P-N)^-PN{P-N). (12) 

The linearization around each of the TLM-constituent pulses (if considered in isolation) carries a pair of zero 
eigenvalues due to the phase invariance. The corresponding cigcnfunction, W, has the shape of the pulse itself (P 
and N, respectively). There is also a generalized eigenfunction [26] for each pair. We expect to have two pairs of 
TLM eigenvalues originating from these pairs of the individual-soliton's ones, the corresponding eigenfunctions being 
(in the lowest approximation) W = P + N and W = P — N. The former one is related to the phase invariance of 
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the system (in other words, to the conservation of the total number of atoms), which means that the corresponding 
eigenvalue pair remains exactly equal to zero, while the other one may become different than zero. Then, it follows 
from Eq. (12) and Eq. (9) that, in the lowest approximation, U = —PN{P — N). In turn, Eq. (11) yields 



K = sign 



r 

,j — ( 



PN{P - Nfdx 



(13) 



Since, by definition, the signs of P(.x) and N{x) are opposite, Eq. (13) implies that the Krein signature of this phase 
eigenmode is negative. A known consequence of this [27] is that an oscillatory instability will be generated by collision 
of this nonzero eigenvalue with other discrete ones (and/or those belonging to the continuous spectrum) that carry a 
positive Krein signature. 

The above consideration indicates that, for every pair of the; interacting 7r-out-of-phase pulses, there exists an 
eigenvalue with negative Krein signature and, consequently, a potentially ensuing Hamiltonian Hopf bifurcation. This 
conclusion fully agrees with our numerical results presented above, that have identified a single oscillatory instability 
for the single TLM soliton, and two such instabilities in the case of bound states of two TLM solitons. 



C. TLM solitons in the absence of the OL 



Finally, it is possible to present analytical arguments justifying the existence of TLM solitons (both single and 
multiple ones) even when the harmonic trap is overwhelmingly stronger than the OL potential (see the numerical 
results shown in the last panels in Figs. 2 and 4), or when the OL potential is simply absent. Within the framework 
of the same perturbation theory which was employed above to derive Eq. (6), each of the pulses constituting the 
TLM state is then subject to the action of two forces. One of them is exerted by the magnetic trap, while the other 
comes from the interaction between the pulses. Accordingly, the effective potential is 



Kff(a, 6; A?!>) = -Srf cos (A^!>) exp 



V2 



(14) 



Predictions following from this approximation were tested against direct numerical results obtained in the absence 
of the OL, i.e., with Vq = in Eq. (1). In Fig. 7, an example is displayed for e = 0.2 and ji = —1.5. The top panel 
shows the effective potential (14), which, in this case, has a minimum at S,i,2 = ±1.35. The middle and bottom panels 
present a numerical solution, with maxima of \u{x)\ at ^i_2 ~ ±1.4. Thus, the perturbative approximation not only 
explains the existence of the TLM state in attractive BECs even in the absence of OL, but also accurately predicts 
the location of centers of the pulses whose concatenation generates the TLM structure. This, in turn, supports the 
numerical findings presented in section II, which suggested that TLM solitons in the attractive BEG would indeed 
persist even for very tight magnetic traps. Thus, we conclude that the TLM solitons are a very robust feature of the 
system. 
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FIG. 7. The top panel shows the eflfective potential V^fi, as a function of the coordinate ^ of the right pulse, which is induced 
by the repulsion between the pulses and the restoring force from the magnetic trap, following the analytical expression (14). The 
potential minimum (equilibrium position) is located at ^ « 1.35. The middle and bottom panels show the TLM configuration 
and its stability, as found in the numerical form. The equilibrium position corresponding to the numerical solution is inferred 
from the middle plot to be at ^ « 1.4. The parameters are Vo = 0, e = 0.2, and /i = —1.5. 

IV. CONCLUSIONS 

In this work we have demonstrated that counterparts of twisted localized modes (TLMs) , which were earlier found in 
the discrete one-dimensional NLS equation, can exist as robust objects in attractive, effectively one-dimensional Bose- 
Einstein condensates (BECs), provided that the condensate is confined by the optical lattice (OL) and/or parabolic 
magnetic trap. The TLMs in BECs may be considered as bound states of two fundamental solitons with a phase 
difference tt between them, or, alternatively, as a dark soliton embedded in a broader bright one. The family of the 
TLM solitons and their stability were investigated in detail, by means of numerical continuation in the strengths of 
the OL or magnetic-trap parameters, as well as in the chemical potential. Stable bound states of two TLMs, alias a 
bound state of three fundamental solitons with a tt phase shift between adjacent ones, were also found. In the cases 
when the TLMs are unstable, the development of the instability was monitored by means of direct simulations, with 
a conclusion that the TLM rearranges itself into a fundamental soliton. 

We have also developed an analytical approach, which treats each fundamental soliton as a quasi-particle subject 
to the action of effective forces exerted by the OL, magnetic trap, and interaction with another soliton. The analysis 
predicts the existence of the stable bound state of two fundamental solitons with a phase shift of tt between them. 
The smallest possible distance between the bound solitons, predicted by this approach, is in good agreement with the 
numerically found size of the corresponding TLM. The analysis can be readily extended to predict that multi-pulse 
bound states with the phase shift tt between adjacent solitons also exist and may be stable. 

The origin of the instability of the TLMs, in the case when they are unstable, was also examined by means of an 
analytical approximation, and was found to stem from Hamiltonian Hopf bifurcations, which occur when eigenvalues 
with negative Krein signature collide with positive-signature ones. This result was based on the existence of a 
negative-Krein-sign eigenvalue set, that was constructed by means of the perturbative analysis. 

It will be very interesting to generalize these considerations to higher dimensions - in particular to two-dimensional 
settings. In that case, a relevant issue is the possible existence of a 2D counterpart of the TLM in the form of a bright 
vortex (which was found in the 2D discrete NLS equation, where it was shown to be stable under certain conditions 
[29]). A very recent preliminary result is that bright vortices, stabilized by a 2D optical lattice, can indeed be found 
[30]. 
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